function [Mom, dataSim, f, entrant, lambda_e, equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon] = moments(bParams,oParams,setts)
% Calculate model's moments

% Variables:
% bParams,oParams,setts		see readme

%% Settings

% Set parameters
	pbar = oParams.pbar;
	n = oParams.n;
    n_e = oParams.n_e;
    product_distribution_entrants = oParams.product_distribution_entrants;
    cdf_product_distribution_entrants = cumsum([0 product_distribution_entrants]);
    theta = bParams.theta;
	m_sim = setts.m;
	T = setts.T;
	T0 = setts.T0;
	dt = setts.dt;

%% Solve Model

% Calcualte update matrix
	upMat = updatematrix(bParams,oParams,setts);

% Calculate equilibrium
    [f, entrant, lambda_e, equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon, C] = equilibrium(bParams,oParams,setts,upMat);
	Nc = length(C);

%% Approximation entry distribution (using steady state distrubtion from the no debt model)
    [approx_steady_state_dist,c_vec] = steadystate_distribution(bParams,oParams,f,C,lambda_mat,equity_mat,optimal_coupon,n_e,product_distribution_entrants,coupon_mat);
    cdf_approx_steady_state_dist = cumsum(approx_steady_state_dist);

%% Simulations

% Initial settings simulations
    rng(13123);
	draws = rand(round(T/dt),m_sim);
    draws_initial_products = rand(round(T/dt),m_sim);
    draws_new_products = binornd(n,theta,round(T/dt),m_sim);
	simulations = zeros(8,round(T/dt),m_sim); %({P0, age, P_t, C_t, default, exit}, Time, # of simulations)

% Loop over simulations
    for i = 1:m_sim
        
% Loop over time
        for t = 1:(T/dt)
            
% Initiation at t = 1
            if t == 1
                p0_dum = sum(draws_initial_products(t,i) >= cdf_approx_steady_state_dist);
                simulations(1,t,i) = mod(p0_dum,pbar+1)-1; % P0
                if simulations(1,t,i)<0 % If the draw = 1 then wrong value for P0 will be assigned, which we need to correct
                    simulations(1,t,i) = pbar;
                end
                simulations(2,t,i) = 0; % Age is zero
                simulations(3,t,i) = simulations(1,t,i); % Pt
                simulations(4,t,i) = c_vec(ceil(p0_dum/(pbar+1))); % Ct
                simulations(5,t,i) = find(C == simulations(4,t,i)); % Ct index
                simulations(6,t,i) = 0; % Default
                simulations(7,t,i) = 0; % Exit
                simulations(8,t,i) = 0; % Refinancing

            end

% Update as if nothing happened (no default and exit as well)
            simulations(:,t+1,i) = simulations(:,t,i);
            simulations(2,t+1,i) = simulations(2,t+1,i) + dt;
            simulations(6,t+1,i) = 0; % Default
            simulations(7,t+1,i) = 0; % Exit
            simulations(8,t+1,i) = 0; % Refinancing

% Transition t -> t+1  
            change = ((draws(t,i) <= lambda_mat(simulations(5,t,i),simulations(3,t,i)+1) * dt + simulations(3,t,i) * f * dt));
            if change == 1
% Update Pt
% Product creation
                if (draws(t,i) > simulations(3,t,i) * f * dt)
                    simulations(3,t+1,i) = min(simulations(3,t,i) + draws_new_products(t+1,i),pbar);
% Creative destruction
                else
                    simulations(3,t+1,i) = simulations(3,t,i) - 1;
                end

% Default and Exit
                if equity_mat(simulations(5,t+1,i),simulations(3,t+1,i)+1) == 0
% Exit
                    if simulations(3,t+1,i) == 0
                        simulations(3,t+1,i) = sum(draws_initial_products(t+1,i) >= cdf_product_distribution_entrants);
                        simulations(7,t+1,i) = 1;
                    end
% Default and creat new firm
                    simulations(1,t+1,i) = simulations(3,t+1,i); % P0
                    simulations(2,t+1,i) = 0; % Age is zero
                    simulations(4,t+1,i) = optimal_coupon(simulations(3,t+1,i)+1); % Ct
                    simulations(5,t+1,i) = find(C == simulations(4,t+1,i)); % Ct index
                    simulations(6,t+1,i) = 1; % Default
                end
            end

        end
    end

%% Generate moments

% Rescale data
	p0sim_g = squeeze(simulations(1,:,:));
	ptsim_g = squeeze(simulations(3,:,:));
	ctsim_g = squeeze(simulations(4,:,:));
	default_g = squeeze(simulations(6,:,:));
    exit_g = squeeze(simulations(7,:,:));
    
	unique_C = unique(ctsim_g); 
	ctsim_g_id2=zeros(size(ctsim_g));

	if setts.nodebt==1 
		ctsim_g_id2 = ones(size(p0sim_g));
	else
		for k=1:length(unique_C)
			ctsim_g_id2(ctsim_g==unique_C(k)) = find(C==unique_C(k));
		end
	end

% Calculate incidence and value 

% Only for non-defaulted firms
	incidence_sim_g = (ptsim_g(2:end,:)-ptsim_g(1:end-1,:))>0;       
	incidence_sim_g = [zeros(1,size(incidence_sim_g,2)); incidence_sim_g];

% Trash all intermediate data
	idxSim = round((1:1/dt:(T+dt-T0)/dt));

	p0sim = p0sim_g(idxSim(2:end)',:);
	ptsim = ptsim_g(idxSim(2:end)',:);
	ctsim = ctsim_g(idxSim(2:end)',:);
	ctsim_id2 = ctsim_g_id2(idxSim(2:end)',:);

	NFirm = size(p0sim,2);
	n_T = size(p0sim,1);

	incidence_sim = nan(length(idxSim)-1,NFirm);
	default_sim = incidence_sim;

	for i=1:length(idxSim)-1
		default_sim(i,:) = sum(default_g(idxSim(i):idxSim(i+1),:));
		incidence_sim(i,:) = sum(incidence_sim_g(idxSim(i):idxSim(i+1),:));
	end

% Other variables
	lambda_sim = lambda_mat(sub2ind(size(lambda_mat),ctsim_id2,ptsim+1));
	equity_sim = equity_mat(sub2ind(size(equity_mat),ctsim_id2,ptsim+1));
	debt_sim = debt_mat(sub2ind(size(debt_mat),ctsim_id2,ptsim+1));
    leverage_sim = (debt_sim./(equity_sim+debt_sim+10^-5));

% Calculate all the moments
	
    mean_lambda = mean(mean(lambda_sim));
    
	exit_rate = sum(sum(default_sim>0))/numel(default_sim);
	
	mean_incidence = mean(mean(incidence_sim>0));
	
	mean_leverage = mean(mean(leverage_sim));
	mean_leverage_f = mean(leverage_sim);
	leverage_sim_w = leverage_sim - repmat(mean_leverage_f,[n_T 1]);   
	var_leverage = var(reshape(leverage_sim_w,1,NFirm*n_T));

	dataSim.pt = ptsim;
    dataSim.pzero = p0sim;
	dataSim.ct = ctsim;
	dataSim.leverage = leverage_sim;
	dataSim.default = default_sim;
	dataSim.incidence = incidence_sim;
	dataSim.debt = debt_sim;
    dataSim.equity = equity_sim;
	dataSim.C = C;
	dataSim.lambda = lambda_sim;

%% Convergence to steady state
    disp("% difference mean firm size beginning and end of simulations: " +  (mean(ptsim(end,:)) - mean(ptsim(1,:)))/mean(ptsim(1,:)) + " and in raw difference " +  (mean(ptsim(end,:)) - mean(ptsim(1,:))))

%% Masses
    M_incumbents = 1/mean(mean(ptsim_g));
    exit_rate_with_zero_products = mean(mean(exit_g))/dt;
    investment_rate_incumbents = mean_lambda * theta * n;
    fi = M_incumbents * investment_rate_incumbents;
    M_entrants = M_incumbents * exit_rate_with_zero_products / lambda_e;
    investment_rate_entrants = lambda_e * (1:pbar) * product_distribution_entrants';
    fe = M_entrants * investment_rate_entrants;
    disp("Approximate difference between rate of creative destruction and number of products created: " +  (M_incumbents * investment_rate_incumbents + M_entrants * investment_rate_entrants - f) )
    flow_new_entrants = M_entrants * lambda_e;

%% Moments table
    Mom = table(f, fe, fi, mean_leverage, var_leverage, mean_incidence, exit_rate, M_incumbents, M_entrants, flow_new_entrants);

end